Non-Nearest-Neighbor Interactions in Nonlinear Dynamical Lattices 
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We revisit the theme of non-nearest-neighbor interactions in nonlinear dynamical lattices, in 
the prototypical setting of the discrete nonlinear Schrodinger equation. Our approach offers a 
systematic way of analyzing the existence and stability of solutions of the system near the so-called 
anti-continuum limit of zero coupling. This affords us a number of analytical insights such as the 
fact that, for instance, next-nearest-neighbor interactions allow for solutions with nontrivial phase 
structure in infinite one-dimensional lattices; in the case of purely nearest-neighbor interactions, 
such phase structure is disallowed. On the other hand, such non-nearest-neighbor interactions can 
critically affect the stability of unstable structures, such as topological charge S = 2 discrete vortices. 
These analytical predictions are corroborated by numerical bifurcation and stability computations. 



INTRODUCTION 

Undoubtedly, the discrete nonlinear Schrodinger 
(DNLS) equation is one of the most fundamental non- 
linear lattice dynamical models. Part of its relevance is 
due to its being one of the main discrete analogs of the 
famous (and integrable in 1+1 dimensions) continuum 
nonlinear Schrodinger (NLS) equation 0,0). The latter 
is encountered in a wide spectrum of applications; it is 
the relevant dispersive envelope wave model for the elec- 
tric field in optical fibers [H, [j] , for the self-focusing and 
collapse of Langmuir waves in plasma physics 0, Q or for 
the description of freak waves (the so-called rogue waves) 
in the ocean [7[. 

On the other hand, the DNLS model is of physical 
interest in its own right within a diverse host of applica- 
tions; see the reviews [H, Q , as well as the recent book 
[Io| . One of the principal directions that spurred a con- 
siderable interest in the DNLS model was the area of the 
nonlinear optics of fabricated AlGaAs waveguide arrays 
llj . In this seting, a plethora of important phenomena 
have been experimentally investigated including discrete 
diffraction, Peierls barriers (the energetic barrier that a 
wave needs to overcome to move over a lattice), diffrac- 
tion management (the periodic alternation of the diffrac- 
tion coefficient) [l2j and gap solitons (structures localized 
due to nonlinearity in the gap of the underlying linear 
spectru m) 1 1 31) among others [14j , More recently, quasi- 
periodic [lo| and even completely disordered [T(| variants 
of such waveguide lattices have been constructed, allow- 
ing the observation of localization type phenomena (even 
within the linear regime for ones of the Anderson type). 
Many of these results and the corresponding theoretical 
activity that they triggered has been summarized in a 
number of reviews see e.g. fl7l [l8j]. 

Another independent and considerably different phys- 
ical setting where DNLS-type models naturally arise is 
that of Bose-Einstein condensates (BECs) trapped within 
periodic (so-called, optical lattice) potentials created by 
counter-propagating laser beams in one-, two- or even 
three- spatial dimensions [3]. This area of atomic 



physics has also experienced a huge growth over the past 
few years, including the prediction and manifestation of 
modulational instabilities (i.e., the instability of spat ially 
uniform states towards spatially modulated ones) [201 ] . 
the observation of gap solitons 2l|, Landau-Zener tun- 
neling (tunneling between different bands of the periodic 
potential) |22j and Bloch oscillations (for matter waves 
subject to combined periodic and linear potentials) [23I ] 
among many other salient features; reviews of the theo- 
retical and experimental findings in this area have also 
recently appeared in [24], [25| . 

Our aim in the present communication is to consider 
a somewhat non-standard variant of the DNLS equation, 
namely one in which the interaction kernel within the lin- 
ear term allows interactions among any set of neighbors. 
Such, so-called, nonlocal versions of the DNLS model 
have been abundantly considered in the past, yielding nu- 
merous interesting features. In particular, it was shown, 
for instance, that, if the interaction strength decays suf- 
ficiently slowly as a function of distance, it can give rise 
to bistability of fundamental solitary waves (centered on 
a single site) [26], which may find applications in their 
controllable switching 27[ ■ Dynamical lattices with long- 
range interactions also serve as models for energy and 
charge transport in biological molecules [28(. Systems 
including the next-nearest-neighbor (NNN) interaction 
have also been proposed as models of polymers (29[, as 
well as in the case of optical waveguide arrays near the 
zero-dispersion point [301 ] . Nonlinear lattice models with 
competing short- and long-range interactions were stud- 
ied not only in one, but also in two dimensions 31 1. 
Lastly, it is relevant to mention that quantum DNLS 
models with nonlocal interactions were another subject 
of study (by means of the Bethe ansatz) [321] . It should be 
pointed out that such longer-range interactions are often 
directly relevant to physical applications. For example, in 
the case of the waveguide arrays, the NNN interaction is 
also, in principle, present yet its relative strength depends 
on the ratio of the separation between adjacent cores to 
the wavelength of light. A specific, zigzag-shaped, model 
of an optical array involving a NNN coupling was in- 
troduced in Ref. [33j , allowing for essentially arbitrary 
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ratios of nearest- to next-nearest-neighbor interactions, 
depending on the specific geometry. Lately, long-range 
interactions, especially of a dipolar type have become 
of considerable interest in BECs due to the significant 
component of dipole-dipole interactions in 52 Cr. These 
effects, as well as the significant component of studies of 
such dipolar BECs in optical lattices have recently been 
summarized in [34(. 

The scope of the present communication is to explore 
the non-local variant of DNLS in the vicinity of the so- 
called anti-continuum (AC) limit, where all individual 
sites (i.e., waveguides in optics or wells of the optical lat- 
tice in BEC) are uncoupled. From there, in the spirit of 
[3^ . we can develop a perturbative expansion using the 
strength of the interactions as the relevant small param- 
eter. This allows us to derive conditions for the persis- 
tence of states near the AC-limit, as well as to formulate 
a framework to address their linear stability. Both the 
existence conditions, as well as the stability conditions al- 
low us to infer valuable conclusions. The existence condi- 
tions suggest (and numerical computations corroborate) 
the possibility that solutions with a nontrivial phase dis- 
tribution may arise in the Id setting, a feature which 
is absent in the case of nearest-neighbor (NN) interac- 
tions; even the sole inclusion of NNN interactions allows 
to achieve this feature. On the other hand, concerning 
stability, the inclusion of NNN interactions may also play 
an important role, as we illustrate with a two-dimensional 
example: the stabilization of a topological charge 5 = 2 
vortex (which is unstable in the presence of NN interac- 
tions) is numerically observed and analytically justified. 

Our analytical considerations are presented in section 
II, while section III is devoted to the corresponding nu- 
merical results (and their comparison to theory). Finally, 
section IV summarizes our findings and suggests a num- 
ber of interesting future directions. 



ANALYTICAL CONSIDERATIONS 

Our nonlocal variant of the DNLS equation will be of 
the following form: 
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where the overdot denotes temporal (in the case of BECs 
in optical lattices) or spatial (in the case of waveguide 
arrays) derivatives, n denotes the site index (waveguide 
index in the optical case, or optical lattice well in BECs), 
and e controls the strength of the coupling among sites, 
which is also determined by the "coupling matrix" with 
elements k nm . The dependent variable u represents the 
complex envelope of the electric field in optics, or the 
complex wavefunction of the BEC in atomic physics. The 



AC limit is defined by e — 0, in which case the unper- 
turbed energy of the uncoupled oscillators reads 



(2) 



Once the inter-site coupling is turned on, the relevant 
perturbation in the energy reads: 
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Notice that at the AC limit, the solutions for the excited 
oscillators will be u n = e l9 ", where the 9 n 's are arbitrary 
phase parameters {u n = corresponds to the non-excited 
sites). 

To determine the persistence of the waves, as was il- 
lustrated in [3^| (see also the general framework of [37]] 
or the nearest-neighbor case of 35|) one has to evaluate 
the perturbed energy at the unperturbed limit solution, 
in which case, we can straightforwardly evaluate it to be: 
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where N is the number of excited (adjacent in this case, 
for simplicity, -although the theory can be appropri- 
ately generalized for non-adjacent) oscillators. The gen - 
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eral persistence conditions then read de n E\ = 
i.e., the unperturbed wave needs to correspond to an ex- 
tremum of the perturbation energy in order to persist 
This leads to the condition 
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On the other hand, the general stability theory of [37 1 
(again see [351 ] for a particular application in the case of 
nearest-neighbor interactions) allows us to quantify the 
bifurcation of the N—l eigenvalues from A = 0. There are 
N eigenvalues at A = when e = 0, and only one of them 
remains there for finite e, due to the U(l) symmetry and 
associated phase invariance, while N — 1 acquire small 
(0(v / e) or smaller) values, as the sites become coupled. 
These eigenvalues satisfy the reduced eigenvalue problem 
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(7) 

(8) 



i.e., the stability is effectively determined by the Jacobian 
of the persistence conditions. It is worthwhile to note 
that this general formulation encompasses the nearest- 
neighbor one of [35j as a special case. 
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The above formulas, namely Eq. ([5]) describing the 
existence and Eqs. ((6])- (5} quantifying the stability un- 
der weak coupling are the main analytical results of the 
present study, which we now proceed to utilize in specific 
applications where additional (to the nearest neighbor) 
interactions exist. 



NUMERICAL RESULTS 

In this section, we apply the above analysis to a few 
select cases where interactions other than purely near- 
est neighbor ones are included, so as to provide a sense 
of the kind of additional features that such settings may 
entail and how these features are brought forth through 
the above analysis. We start with what arguably can 
be dubbed the simplest extension of NN interactions, 
namely the inclusion of the NNN ones which is directly 
physically realizable, as per the setting of [33l ]. 

We start from the case of equal NN and NNN inter- 
actions only, i.e., k nm = 5 m ,„±i + S m . n ± 2 , where 6 de- 
notes the Kronecker 5. Considering then a N = 3 site 
configuration (the smallest one that "perceives" of the 
presence of the NNN interaction), it is straightforward 
to intuitively realize that this is essentially the analog 
of an equilateral triangle where all pairwise interactions 
exist and are equal to each other. This is a situation 
that has been considered both in DNLS [3£l (see also the 
early work of and in Klein-Gordon [40, E| type set- 
tings, motivated by applications of triangular lattices in 
photorefractive crystals and dusty plasmas. It has been 
illustrated therein that the fundamental linearly stable 3- 
excited-site configuration is that of a so-called vortex of 
topological charge S = 1, namely a situation where the 
phase runs from at the first excited site, to 2tt/3 at the 
second one and 4ir/3 at the third one. Remarkably, all 
configurations with the standard (from the NN setting of 
[IB)]) case of or -k phases, are unstable. Equally impor- 
tantly the NN case, as per the arguments of [35(, never 
allows for the existence of phases other than or n in 
Id configurations. Hence, one of the principal features 
of these non-nearest-neighbor interactions is that they 
may enable configurations that would not be possible in 
standard NN settings. Moreover, as indicated above in 
this case, the eigenvalues can be analytically computed 
(through a discrete Fourier transform applied to the "tri- 
angle" of sites (which has periodic boundary conditions). 
This can be found to lead to the eigenvalues 



Xj = ±W8ecos(A(9)sin 2 ( ^ 
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for j = 1, 2, 3. The stable case features the relative phase 
between adjacent sites AO = 2ir/3 and yields a double 
pair of imaginary eigenvalues ±\/3ei which are compared 
to the full numerical results in Fig. [1] (also the analytical 
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FIG. 1: The figure shows a 3-site solution with a nontriv- 
ial phase distribution (0, 2-k/S, A-k/3) in the Id case with 
equal nearest and next-nearest-neighbor interaction strengths 
of unity, i.e., k n , n ±2 = k n , n ±i = 1- The first two panels in the 
top row show the real and imaginary part of the solution for 
e = 0.15, while the bottom left shows the corresponding spec- 
tral plane (A r , Ai) of the linearization eigenvalues A = X r +iXt. 
The second panel of the bottom row shows the relative phases 
A9 = 62 — 9i (blue solid line) and #3 — 81 (red dashed line), 
which does not change significantly from its anti-continuum 
limit as e is increased. Lastly, the top right and bottom right- 
most panels show the maximal real eigenvalue (larger than 
zero -and hence giving rise to an instability- for e > 0.116) and 
lowest imaginary eigenvalues, respectively (blue solid lines). 
The theoretical prediction for the latter is given by the red 
dashed line. 



existence prediction for the relative angle among adja- 
cent sites is tested through the numerical results of the 
figure). Good agreement is obtained between the ana- 
lytical and numerical results in both existence and sta- 
bility, for small values of e, while, as expected, relevant 
deviations increase for larger values of the perturbation 
parameter. We also note in passing that the non-trivial 
phase vortex structure becomes unstable eventually in an 
oscillatory way due to a complex quartet of eigenvalues 
for e > 0.116 as also illustrated in Fig. [TJ due to the col- 
lision of one of the imaginary pairs with the continuous 
spectrum extending above X = i (and below A = — i). 

To illustrate the generality of our analytical consid- 
erations in settings other than the special, degenerate 
case of Fig. [1] we considered another physically realiz- 
able (as per [33j) case where k n n ±2 — 0.7fc n-n= |-i = 0.7 
(while all other interactions are absent). Again, focus- 
ing on the 3-site interaction, we realize that the three 
persistence conditions can be explicitly solved, yielding 
two possible scenarios. Either sin(#„ — 6 m ) = for all 
combinations of n, m G {1, 2, 3} or B\ — 62 = #2 — $3 and 
cos(#i — 8 2 ) = —k n . n+ i/(2k ntn+ 2)- Notice that this en- 
compasses the discussion of equal strengths of NN and 
NNN interactions as a special case and additionally re- 
veals the condition under which the nontrivial phase 
distributions discussed above can generally materialize, 
since it is necessary that \k n ^ n+ i/(2k n . n+2 )\ < 1- In the 
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FIG. 2: An example of weaker (than nearest-neighbor) next- 
nearest-neighbor interaction strength with fc n ,n±2 = 0.7 while 
k n ,n±i — 1: the phase distribution is shown on the left (and its 
anti-continuum limit, as theoretically predicted is depicted by 
squares), while the predictions of the linear stability analysis 
(red dashed lines) are compared with the full numerical results 
(blue solid lines) as a function of e on the right. The real 
part of the corresponding eigenvalues becomes nonzero for 
e > 0.109 (as shown in the top right panel). 



case of Fig. [3J the analytical prediction yields a phase 
difference of 0\ — O2 = 2.3664 {9\ — 83 is twice that) in 
excellent agreement with the numerical observations. No- 
tice also that the relative phases appear to depend only 
weakly on e in the continuations over the latter parame- 
ter shown in Fig. O Additionally, our stability considera- 
tions allow us in this case to obtain the analytical predic- 
tions for the two (now split due to the asymmetry in the 
interactions) pairs of small eigenvalues A = ±v4~2858ei 
and A = ±y/1.37Uei. Good agreement is once again 
found with the full numerical linear stability results for 
small couplings e, while this agreement deteriorates when 
increasing the coupling. The solution becomes unsta- 
ble once again through the oscillatory instability arising 
upon collision (for e > 0.109) of the larger one among 
the above 2 pairs with the continuous spectrum. 

However, the above features such as the nontrivial 
phase distribution are not a unique feature of cases with 
just NN and NNN interactions. As an alternative ex- 
ample that shares these features, we consider the case 
where k nm = 5 m . n ±i + <S m , n ±2 + <S m ,n±3- Contrary to 
what might be naively expected, we do not find here the 
generalization of a vortex including 4-sites. While the 
lowest order persistence conditions are satisfied by such 
a solution, higher order ones are not. However, this does 
not preclude the existence of genuinely complex solutions 
with a complicated phase structure. An example of this 
kind involving 7 sites is shown in Fig. [3] The phases 
of the 7 excited sites are found at the AC limit to be 
#1=0 (since one of them can always be chosen arbitrar- 
ily), 9 2 = 3.2831, 3 = 1.8780, 4 = 4.0894, 9 5 = 0.01766 
and 07 — 1.8957. Then, the corresponding eigenvalue 
predictions from the reduced eigenvalue problem are: 
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FIG. 3: A more complicated Id example of a 7-site wave with 
nontrivial phase distribution for k n ,n±i = k n ,n±2 = fcn,n±3 = 
1. The top left panels show the real and imaginary parts 
of the solution and its spectral plane for e — 0.2; the top 
right shows the relative phases A6>; = 9, — 8\ for i — 2, ... ,7 
(i = 2: blue solid, i = 3: red dashed, i = 4: green dash- 
dotted, i = 5 thick blue solid, i — 6: thick red dashed, i — 7: 
thick green dash-dotted); the theoretical anti-continuum limit 
is given by squares. Lastly, the bottom panels show the six 
relevant eigenvalues (blue solid lines) and the corresponding 
theoretical predictions (red dashed lines). The real (top) and 
imaginary (bottom) parts of the eigenvalues are given - it is 
clear that the solution becomes unstable for e > 0.062. This 
inset contains a blowup of the imaginary parts for small e to 
showcase the agreement between numerics and theory. 



A = ± V7.0792e i, A = ± V6.875d , A = ±\/ 4.2006e i, 
A = ±V3.4994d, A = ±V0.2248ei, A = ±V0.1334ei. 
These are compared to the full numerical eigenvalues 
of the problem (as are the existence results) in Fig. [3] 
Despite the complexity of the eigenvalue structure, the 
quality of the agreement between theory and numerics is 
clear from the eigenvalue inset and the phase comparisons 
(near the AC limit). Notice that the solution becomes 
unstable for e > 0.062 due to collision of the largest of 
the above pairs with the continuous spectrum and in- 
creasingly so due to additional collisions for e > 0.088, 
e > 0.121 and e > 0.128. 

As our last example and in order to demonstrate a 
case where the non-nearest-neighbor interactions have a 
direct impact on stability, we consider a prototypical two- 
dimensional case of an unstable structure with vorticity, 
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FIG. 4: An example of the stabilized S — 2 vortex for e = 0.02 
(the real contour of the solutions is shown in the top left and 
the imaginary one in the top right). The spectral plane in- 
dicating the stability of the solution is shown in the bottom 
left while the corresponding numerically obtained (imaginary) 
eigenvalues are in the bottom right (blue solid, as compared 
with the red dashed theoretical prediction). The solution be- 
comes unstable for e > 0.032. 



namely the vortex of topological charge S = 2. It is 
well-known 3^, 4^] that this structure is unstable and a 
considerable effort has been invested in stabilizing it ei- 
ther via defects |43l| or geometrically through a rhombic 
configuration [4J]. Here, we illustrate that NNN inter- 
actions can also play a stabilizing role for this coherent 
structure. In particular, we examined a special case (but 
verified that the result was also robust for other choices) , 
where k n m = 1 for |n— m| = 1, k n m — 1 for |n— m| = 2, 
while fc nm — 2 for |n — m| = \/2 (where the last choice 
was made because of the two lattice paths connecting the 
relevant neighbors). In this case, (but also with other 
choices of the relative weights), it is straightforward to 
show that the square configuration with 8j — jn/2, where 
j = 1, . . . , 8 remains a solution and that the correspond- 
ing reduced eigenvalue problem has a single eigenvalue 
pair of A = ±Ay/ei, another single pair of A = ±Vfei, 
and a quadruple pair of A = ±2y / ei (while one pair is ex- 
actly at and one is to leading order but also remains 
imaginary at higher orders). This stabilized variant of 
the S = 2 vortex can be observed in Fig. QJ where again 
the analytical predictions are compared to the numerical 
results, yielding increasingly good agreement once again 
for smaller values of e. Notice that the solution becomes 
destabilized for e > 0.032, while additional quartets arise 
for e > 0.048, e > 0.074 and e > 0.087. 



CONCLUSIONS 

In the present work, we have revisited the theme of 
non-nearest-neighbor interactions in settings of the dis- 
crete nonlinear Schrodinger type motivated by applica- 
tions especially in optical (but also in other areas such 



as atomic) physics. We illustrated that an analytical ap- 
proach enables a number of interesting conclusions and 
possibilities that did not seem to have been realized be- 
fore, to the best of our knowledge, both as concerns the 
existence problem on the infinite lattice, and as concerns 
the stability properties in the context of localized solitary 
wave solutions. 

In particular, we focused on the vicinity of the anti- 
continuum limit of vanishing coupling. There, we could 
obtain explicit conditions for the persistence of localized 
modes, as well as analytical predictions for the eigenval- 
ues of linearization around such solutions. This analy- 
sis allows to conclude that the non-nearest-neighbor set- 
ting has a number of special properties. More specif- 
ically, even just next-nearest-neighbor interactions (but 
also more general ones) create the possibility for solutions 
with more complicated phase profiles than just phases of 
and 7r as was the case for nearest-neighbor interactions 
[3Fjj ]. On the other hand, these longer range interactions 
can play a critical role in the stability of solutions that 
were found to be unstable in the nearest-neighbor set- 
ting. We presented a particular such example in the case 
of the topological charge S = 2 discrete two-dimensional 
vortex. 

Although our work lays the ground for understanding 
some of the main features of these non-nearest-neighbor 
settings, a number of interesting questions are still out- 
standing. In particular, it would be interesting if some- 
thing could be said more generally about the type of so- 
lutions (and phase distributions) afforded by interactions 
of different range/decay properties. On the other hand, it 
would be interesting if a general classification of stability 
could also be made for different types of intereaction ker- 
nels, in higher or even in one dimension. These questions 
are presently under investigation and will be reported in 
future publications. 
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